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A fast Variational Gaussian Wave-packet method: Size-induced structural transitions 

in large neon clusters. 

lonut Georgescu and Vladimir A. Mandelshtam 
Chemistry Department, University of California, Irvine, CA 92697, USA 

The Variational Gaussian wavepacket (VGW) method is an alternative to Path Integral Monte- 
Carlo (PIMC) for the computation of thermodynamic properties of many-body systems at thermal 
equilibrium. It provides a direct access to the thermal density matrix and is particularly efficient 
for Monte-Carlo approaches, as for an A/'-body system it operates in a non-inflated 3N dimensional 
configuration space. Here we greatly accelerate the VGW method by retaining only the relevant 
short-range correlations in the (otherwise full) 3N x 3N Gaussian width matrix without sacrificing 
the accuracy of the fully-coupled VGW method. This results in the reduction of the original 0{N^) 
scaling to 0{N^). The Fast- VGW method is then applied to quantum Lennard- Jones clusters with 
sizes up to iV = 6500 atoms. Following Doye and Calvo [JCP 116, 8307 (2002)] we study the 
competition between the icosahedral and decahedral structural motifs in Ncat clusters as a function 
ofN. 
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Introduction 

The Variational Gaussian wavepacket (VGW) 
method [TJ [2] was introduced recently as an alternative 
to Path Integral techniques for the estimation of thermo- 
dynamic and structural properties of large many-body 
systems at thermal equilibrium. The VGW method 
provides a direct and numerically efficient estimate of 
the density matrix e~^^, and particularly its diagonal 
elements p{x) = {x\e~^^\x). As such it has been 
combined successfully with Monte-Carlo techniques 
for computations of thermodynamic and structural 
properties of atomic and molecular clusters [3]-[8]. When 
compared with accurate Path-Integral Monte-Carlo 
(PIMC) calculations for certain neon clusters, the VGW 
method yielded practically identical results [H [9j . The 
Thermal Gaussian Molecular Dynamics [TO",Tr (TGMD) 
was built on top of the VGW for the estimation of time 
correlation functions in quantum many-body systems. 
Other quantum dynamics approaches, such as the Full 
Wigner dynamics [12 , Equilibrium Liouville Dynamics 
[13l [14] also take advantage of the analytically conve- 
nient representation of the density matrix in the VGW 
formalism. 

Given an A/'-body system, the most accurate version 
of the VGW method, the so-called Fully-Coupled VGW 
(FC-VGW) utilizes the full 3N x 3N Gaussian width ma- 
trix G, with the matrix-matrix multiplication step being 
the numerical bottleneck. The 0{N^) numerical scaling 
of the latter limits the applicability of the FC-VGW to 
systems of relatively small size. The so-called Single- 
Particle VGW (SP-VGW) utilizes the block-diagonal 
form of G, each particle represented by a 3 x 3 block. Con- 
sequently, the matrix-matrix multiplication is no longer 
a numerical bottleneck in the SP-VGW, and for systems 
with short-range pair potentials the overall scaling of the 
SP-VGW can be reduced to linear in A" (after appropri- 
ate potential cut-offs are implemented). The drawback of 
the SP-VGW may be an uncontrollable loss of accuracy 
(compared to the FC-VGW). As shown below, only rela- 



tively short-range correlations in the G-matrix are signif- 
icant, even for systems with quite pronounced quantum 
character, such as {para-}l2)N clusters. Here we exploit 
this circumstance and present an improved version of the 
VGW method, called Fast- VGW, which reduces the over- 
all numerical scaling of the FC-VGW from 0{N^) to at 
least 0{N'^) by retaining only the physically relevant off- 
diagonal elements of the (otherwise full) G-matrix, with 
practically no accuracy tradeoff. (Depending on the sys- 
tem, further reduction in the numerical scaling can be 
achieved by applying additional cut-offs and utilizing the 
sparsity of the Hessian matrix.) 

In the next section the VGW formalism is introduced. 
We then analyze the accuracy of the Fast- VGW approach 
by applying it to the ground state calculations of sev- 
eral Lennard- Jones (LJ) clusters for a range of quan- 
tum characters, from moderate (Ne) to relatively strong 
(para-H2). The last section presents an application to 
very large Ncat clusters with up to A" ~ 6500 atoms, 
where the competition between the icosahedral and dec- 
ahedral structural motifs for the ground state is studied 
as a function of the cluster size A". In particular, this 
study provides an improved estimate of the icosahedral- 
decahedral crossover size for Ncat clusters compared to 
that by Calvo and Doye [15 , who used the Harmonic- 
Superposition Method (HSM) to study the structural and 
thermodynamic properties of quantum LJ clusters for a 
range of quantum parameters. 



The Variational Gaussian Wave-packet 
Approximation 

The matrix elements of the density matrix 

p{x,x';f3) := /x\e-^^\x'\ = {x; f3/2\x'; f3/2) (1) 
can be expressed in terms of the wave-packets |x,r) 

\x;T)=e-^^\x), (2) 



which are solutions of the imaginary-time Schrodinger 
equation (or Bloch equation): 



d_ 



\x;t) =H\x;t). 



(3) 



The VGW approach approximates these wavepackets 
with Gaussians 



(r|x;r) « {2n)-^^'^\\G\\-^'^ 



(4) 



X exp 



-Ar-qfG-\r-q)+^ 



where the Gaussian parameters are G = G{r) G R^^^^^, 
the real symmetric and non-negative Gaussian width 
matrix, q = q{T) G M^^, the Gaussian center, and 
7 = l{^) ^ ^7 "the scale factor. 

With this ansatz the diagonal element of the density 
matrix, or just density p(x; /3) := p(x,x;/3), becomes a 
simple expression in terms of G and 7: 



p{x]l3) 



o27(/3/2) 



(4^)3^/2||G(/3/2)||V2 



(5) 



A variational principle [2l |9] provides the equations of 
motion for the Gaussian parameters </(t), G{t) and, re- 
spectively, 7(t) 
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Tr [(VV^f/)G] - (f/) 
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which are propagated with the initial conditions 

q(0) = X, G(0) = 0, 7(0) = 



(6a) 
(6b) 
(6c) 

(7) 



from r = to T = /?/2. By {U), {VU) and (VV^C/) we 
defined, respectively, the expectation values of the total 
energy, its gradient and its Hessian, e.g. 



(f>> 



{x-p/2\U\x-p/2) 
{x;p/2\x;p/2) ' 



(8) 



In the zero temperature limit p 
imizes the functional 



00, the VGW min- 



E 



{x-p/2\H\x-p/2) 
(x;/3/2|x;/3/2) 



dp 



lnp(x; p) 



(9) 



and thus provides a way for estimating both the ground 
state energy and the density. 

The VGW is exact in the harmonic and classical (high 
temperature) limit; the analytical solution of the multi- 
dimensional harmonic oscillator is provided in Ref. pT] . 
Moreover, if U consists of only pair interactions and if 
these can be fitted in terms of Gaussians, then (t/), (V/7) 
and, respectively, (W^U) can be expressed analytically 



and as such are easy to compute. Fortunately, most of 
the pair potentials used in practice, such as the LJ or 
Coulomb, or the Silvera-Goldman [16] potentials can be 
fitted very accurately using a small number of Gaussians. 
The quantum character, better said, the degree of 
quantum delocalization of a system is conveniently de- 
scribed by the de Boer parameter 



A 



h 



(10) 



which relates the quantum delocalization of the wave- 
functions, estimated by the de Broglie wave-length 
hj ^Jme^ to the characteristic length a. The latter charac- 
terizes the range and e, the strength of the pair potential, 
for example the LJ potential is 



f/(ri,) = 4e 



12 



(11) 



The classical limit is obtained for A <C 1. Typical values 
are A = 0.01 for Xe, A = 0.03 for Ar, A = 0.095 for 
Ne and A ~ 0.3 for P-H2. Rewriting the VGW equations 



( 6a 6c ) in reduced units leads to a general set of equations 
with only one free parameter. A: 
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-G{WU) 


(12a) 


dr 


-G(VV^f/)G + A2 


(12b) 
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-^Tr[(VV'^f/)Gl -{[/) 
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(12c) 



This set of equations can be used to study a wide range 
of species by varying just one parameter, A. 



The Fast- VGW method: Numerical tests. 

In the most favorable case depicted so far, namely 
when U consists of pair interactions that can be fitted 
with Gaussians, the computation of the expectation val- 
ues of the energy, its gradient and, respectively, its Hes- 
sian require 0{N'^) operations. Various algorithms, such 
as tree-codes [T7|, can reduce this effort to 0{N\ogN) or 
even 0{N) [18]. With short-range potentials, such as LJ 
or Silvera-Goldman, the 0{N) scaling can be achieved 
easily with a potential cut-off. 

Yet, the dominant operation in Eqs. ( 12a|12c ) is the 
matrix-matrix multiplication in Eq. (12b) which scales as 
0{N^)^ assuming both the Hessian and the G matrix be- 
ing general, dense matrices. This sets a practical limit on 



the size of the system to be treated numerically using the 
VGW method. (Note also that Eq. ^ requires calcula- 
tion of the determinant of G, which also scales as 0{N^) 
for a full matrix, but this calculation is only performed 
once in order to compute the density at r = /3/2.) For 
sufficiently small values of the quantum parameter A, the 
fully-coupled Gaussian wave-packet can be approximated 




FIG. 1: The G matrix for Nei47 at a temperature of T — 
3.56i^ (T = O.le) within the fuh, 3iVx3A^, Gaussian approach 
(FC-VGW). The color density is proportional to the decimal 
logarithm of the matrix elements. 



with product of single-particle Gaussians, thus reducing 
G to a 3 X 3 block diagonal form. Refs. [2] |9] give some 
idea on how the SP-VGW performs for Ne (A ^ 0.1) clus- 
ters, by comparing it to both the FC-VGW and PIMC re- 
sults. While the SP-VGW seems to correctly describe the 
thermodynamics of Ne clusters, it may actually fail, even 
for such weakly quantum systems, to adequately char- 
acterize different cluster configurations with very close 
energy values. For example, this happens for the Nesg 
cluster, for which the FC-VGW is still very accurate [9]. 
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FIG. 2: The ground state energy per particle for a sequence 
of Lennard- Jones clusters as a function of size. The quantum 
parameter is A = 0.10, which approximatively corresponds 
to Ne. Full line: fully-coupled ?>N x 3 A/" Gaussian, dashed 
blue line: the single particle VGW; open circles: the Fast- 
VGW with a correlation distance rcorr = 1.5cr; diamonds: the 
same, but with rcorr = 1.8cr. The shading at the bottom 
indicates the structural motifs of the ground stated : liquid 
(L), Mackay icosahedron (I^), anti-Mackay icosahedron (I^). 

With increasing quantum parameter, particle correla- 
tions gain significance and the off-diagonal blocks of the 



G matrix should not be neglected. We will thus include 
the off-diagonal blocks, but only for pairs of particles, less 
then a certain distance apart, Xij < rcorr, which makes 
G sparse but preserves its positive-definiteness. Such an 
approach is also motivated by Fig. [l] which shows the 
G matrix for a Nei47 cluster at T = 3.b6K (T = O.le) 
within the fully coupled Gaussian (FC-VGW) approach 
(note the logarithmic color scale). Most elements of the 
G-matrix are seven orders of magnitude smaller then the 
dominant ones and can therefore be neglected. 

A first example of the sparse G approach (Fast- VGW) 
is given in Fig. |2J which shows the ground state en- 
ergy per particle Eq/N for a sequence of (LJ)7v clus- 
ters as a function of size. The quantum parameter is 
A = 0.1, corresponding roughly to neon. The ground 
state ener gy was c omputed from Eq. ([9| after propagat- 
ing Eqs. ( 12a|12c ) to (3 = 100 (the Gaussian parame- 
ters become stationary at /3 ~ 50 already). The en- 
ergies obtained with fast- VGW (symbols) are indistin- 
guishable from that using the full 3N x 3A^ G matrix 
(solid black line). Note that the LJ interaction is typ- 
ically cut off at radius Tcutoff = 2.75cr, which is much 
larger than rcorr = l-5cr (open circles). For N — 147, 
the latter value results in a G matrix with 7.2% non-zero 
elements, i.e. the sparsity of G can already make a no- 
ticeable difference. Although the single- Gaussian result 
(the dashed blue line) is displaced, but it runs parallel to 
the FC-VGW result and is able to correctly characterize 
different configurations according to their energies. The 
shaded areas indicate the structural motif of the ground 
state according to the n — A phase diagram reported in 
Ref. [7]. 
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FIG. 3: Same as Fig. |2]but for A 
para-hydrogen. 



0.30, corresponding to 



Fig. [3] shows the same analysis, but for a larger value 
of the quantum parameter, A = 0.3, which corresponds 
to para-H2. The displacement of the SP-VGW result 
becomes irregular. The Fast-VGW using rcorr = l-5cr 
is more accurate, but there is a discrepancy for smaller 
clusters. As indicated by the shading, the ground state 
structures for this size range have liquid-like character. 



meaning that the pair distribution function does not im- 
mediately vanish to zero after the nearest neighbor peak. 
Increasing the correlation cut-off to rcorr = 1.8<j removes 
the discrepancy. The increased correlation cut-off did not 
change the results for the lower A in Fig. [2] 
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FIG. 4: Average time (in arbitra ry units) needed to eval- 
uate the right hand side of Eqs. ( |12a|12c ). Black squares: 
fully-coupled Gaussians; dashed blue: single-particle Gaus- 
sians; open red circles: Fast-VGW with correlation distance 
rcorr = 1.5cr; crosses: Fast-VGW with rcorr = 1.8a. The 
inset compares the running times of the FC-VGW and Fast- 
VGW implementations normalized by the running time of the 
SP-VGW one, underscoring the 0{N^) running time of the 
Fast-VGW implementation. 



Fig. |4] compares the computational cost of the 
three methods presented here: the FC-VGW (full G- 
matrix), SP-VGW (block-diagonal G-matrix) and Fast- 
VGW (sparse G-matrix). The inset shows the running 
times for the FC-VGW and the Fast-VGW implementa- 
tions normalized by the running time of the SP-VGW im- 
plementation. The latter is obviously the fastest, but the 
Fast-VGW is slower by a constant factor and falls in the 
same class of 0{N'^) scaling. No cut-off rad ius w as used 
in estimating the Hessian (VV^/7) in Eq. ( |12b ), which 
is 0{N^) and by far the dominant part of the SP-VGW. 
The linear increase in the inset confirms the 0{N^) cost 
of the full matrix approach. One should stress that 
the Hessian (VV^t/) is still a dense matrix. Applying 
the usual LJ cutoff Tcutoff = 2.75cr or an even larger, 
more "generous" one, would make the Hessian sparse 
too, which would further improve the performance. The 
additional cutoff was not applied here to keep the im- 
plementation as simple as possible. One should also note 
that the increased correlation cut-off did not increase the 
running time significantly. 



Size-induced icosahedral-decahedral transition in 
large neon clusters. 

The Fast-VGW is so far able to provide the same ac- 
curacy as the fully coupled VGW at a fraction of the 
computational cost. Following Ref. [15 we applied 
this new powerful tool to a number of Ncat clusters 
with up to N=6500 atoms to study the competition be- 
tween the structural motifs of the ground state struc- 
ture. The ratio of surface atoms decreases with increas- 
ing cluster size and therefore, a switch from the pre- 
dominantly icosahedral motif of small clusters to motifs 
that optimize the bulk energy over the surface energy 
is expected. Ref. [15] has estimated the crossover sizes 
within the Harmonic-Superposition Method (HSM) for 
both classical and quantum clusters. We expect to im- 
prove the quantum mechanical result of the HSM because 
the VGW accounts inherently for anharmonic contribu- 
tions and relaxes, at the same time, the structure of the 
cluster as the quantum particles become "larger" (more 
delocalized). Since the VGW method is variational and 
is exact for a harmonic potential, it should always pro- 
vide better upper-bound energy estimates than the HSM 
method. 

A detailed structural analysis of large clusters as a 
function of their size is practically unfeasible even in the 
classical case, as the number of local minima increases 
exponentially with size, and is already enormous for sev- 
eral tens of atoms; Global optimization methods have not 
broken the TV = 1000 barrier yet [15l [T9H22] . Moreover, 
except for a better accuracy of the crossover size, a de- 
tailed picture would not contribute much to understand- 
ing the physics of the crossover. Thus, as in Ref. [15 , we 
have also settled for a coarse-grained picture, following 
the sequence of lowest energy configurations of Mackay 
icosahedra and Marks decahedra. The icosahedral motif 
actually favors truncated structures for N > 923, corre- 
sponding to the "magic" number sequence minus the 12 
vertices. Consequently, We have studied N = 911, 1403, 
2045, 2857, 3859, 5071 and 6513. Similarly, for large sizes 
the Marks decahedra favor truncated structures too, so 
we have studied TV = 1103, 1660, 2377, 3274 and 4371. 
We evaluated the energ ies according to Eq. ^ after prop- 
agating Eqs. ( 12a||12c ) to (3 = 100. Within each struc- 
tural type, the obtained ground state energy values are 
then used to interpolate the energy as a function of N 
according to: 



E{N) = asN + bsN^^^ + ceN^^^ + d^ 



(13) 



where qe denotes the bulk and Be the surface contribu- 
tion, respectively. 

Fig. [5] shows the computed ground state energies for 
the icosahedral and the decahedral motifs of (Ne)Ar clus- 
ters estimated by means of the Fast-VGW and of the 
Harmonic-Superposition Method (HSM). The VGW en- 
ergies are systematically higher (HSM is not variational) 
and provide a correct upper bound to the ground-state 
energy. Being exact in the harmonic limit, the VGW 
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The crossover is almost imperceptible on an absolute 
energy scale. The lower panel of Fig. |5] shows the rela- 
tive energy of the decahedral sequence with respect to the 
icosahedral sequence. The relative energy of the icosahe- 
dral sequence to itself is obviously zero. Using the inter- 
polation (13), we estimate the crossover at A^ ~ 3100. 



This is significantly smaller than the quantum HSM re- 
sult N = 4640 and underscores the significance of the 
anharmonic effects. Fig. [5] also shows the energy of 
the classical decahedral motif relative to the classical 
icosahedral motif. The crossover occurs much earlier, 
at A/" ~ 1690. As shown in ref. [l^, the vibrational 
frequencies are lower for icosahedral structures than for 
decahedral ones. With increasing quantum parameter 
A, the zero-point energy of the vibrational modes stabi- 
lizes the icosahedra with respect to decahedra and pushes 
the crossover to larger N. Qualitatively identical behav- 
ior has also been observed with the unoptimized (non- 
truncated) structures (not shown here). The classical 
crossover is N ^ 2200, the VGW predicts N ^ 5400 and 
the HSM estimate is TV - 10000. 
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FIG. 5: The ground-state energy per particle Eq/N of Mackay 
icosahedra versus Marks decahedra estimated by the Fast- 
VGW and the HSM approach. Filled symbols: icosahedra 
(Ih); empty symbols: Marks decahedra (M-Dh). Upper panel: 
absolute energies. Lower panel: the energy of the decahedral 
motifs relative the icosahedral ones; in this panel the lines 
show interpolations according to Eq. 



(13) 



would have yielded the HSM result, had the anharmonic 
contribution been insignificant. One should also note 
that the VGW relaxes the cluster structure, which typi- 
cally increases in size slightly due to the quantum delo- 
calization. The classical result is also included to quan- 
tify the zero-point energy effect. The largest structures 
investigated here are illustrated in Fig. [6J 
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FIG. 6: The four largest structures investigated in this work 
for each of the structural motifs. Upper row: truncated 
Marks-decahedra; lower row: truncated icosahedra. From 
Ref.JJj. 



Conclusions. 



We have introduced a rather straightforward, but prac- 
tically significant improvement of the VGW method 
based on discarding long-range correlations in the rep- 
resentation of the (otherwise fully-coupled) Gaussian 
wavepacket. The computational advantage is tremen- 
dous, reducing the overall cost from 0{N^) to 0{N'^) 
with no visible accuracy tradeoff. The new method al- 
lowed us to study neon clusters with up to 6500 atoms. 
Consequently, we were able to directly investigate the 
competition between the Mackay-icosahedral and Marks- 
decahedral structural motifs when the cluster size is var- 
ied. Additional sparsity in the Hessian could also be ex- 
ploited, reducing the numerical scaling even further, to 
0{N). Our implementation can easily take advantage of 
existent massively parallel sparse matrix packages, and 
as such can be applied to systems with tens of thousands 
of particles. 
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